/*
 *  mecsystem.h
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 15/11/07.
 *
 */

#ifndef _mecsystem_h_
#define _mecsystem_h_

#include "BODIES/body.h"
#include "JOINTS/joint.h"
#include "point.h"
#include "BODIES/flexchain.h"
#include "BODIES/flexfem.h"
#include "BODIES/lattice.h"
#include "BODIES/masslinkchain.h"
#include "BODIES/particlebody.h"
#include "BODIES/rigidbody3d4p.h"
#include "BODIES/ground.h"
#include "BODIES/thermodamper.h"
#include "BODIES/thermostickslip.h"
#include "JOINTS/encastre.h"
#include "JOINTS/encastrepenalty.h"
#include "JOINTS/glue.h"
#include "JOINTS/rigidlink.h"
#include "JOINTS/contact.h"
#include "JOINTS/rigidtetra.h"
#include "JOINTS/torsionfreebeam.h"
#include <iostream>
#include <fstream>
#include <string>
#include <map>
#include "commonutil.h"

using namespace std;

class Process;

class MecSystem
{
public:
	string name;
	string type;
	int dim;
	map<string, Body*, less<string> > mBodies;				/** Map of bodies. Key is name */
	map<string, Joint*, less<string> > mJoints;				/** Map of joints. Key is name */
	map<string, vector<double>, less<string> > pointNames;	/** Initial collection of points. Existence is not checked. */
	map<string, Point, less<string> > mPoints;				/** Map of points belonging to the system */
	vector<string> childSysNamesG;							/** Names of systems within the system (global) */
	map<string,MecSystem*,less<string> > mChildSys;			/** Pointers to systems within the system, local to the process */
	string parentSysName;									/** Parent system name*/
	Process* parentProc;									/** Parent process */
	bool isWellDefined;
	map<string,string,less<string> > mChangedPoints;
	map<string,vector<string>,less<string> > pointSets;
	map<string,vector<double>, less<string> >	mPointIniPos;	/** Points with position initial conditions*/
	map<string,vector<double>, less<string> >	mPointIniVel;	/** Points with velocity initial conditions*/
//	map<string,vector<double>, less<string> >	mJointIniPos;	/** Joints with position initial conditions*/
			
public:	
	MecSystem();
	MecSystem(string, ifstream &, int, Process*, ofstream &);			/** Standard constructor */
	MecSystem(string, map<string,MecSystem*,less<string> >, Process*, int);		/** From a collection of Systems */ 
	~MecSystem()
	{}
	
	void readStandard(ifstream &, int);				/** read standard type system */
	void checkCreatedSystem();						/** checks that it has at least one body or one system.  Formal check only */
	MecSystem* getPTSystem(string);					/** get a pointer to any system, belonging or not to the same process */
	MecSystem* getPTSystemLocal(string);			/** get a pointer to any system of the local process */
	Point* getPTPoint(string);						/** get a pointer to any point, belonging or not to the same process, based ONLY on the name */
	Point* getPTPointLocal(string);					/** get a pointer to any point local to the process */
	void setParenthood();							/** Parenthood relations within the local process */
	void createPoint(const Point*);					/** Creates and adds one point to the system */
	void createPointAtRoot(const Point*);			/** Creates and adds one point to the root system */
	void createPoints();							/** Actually creates the points, once systems, bodies and components are defined */
	void addPointName(const string&, const vector<double>&);	/** Add point names to the map pointNames */
	void setPPointersCompBodies();						/** Set the pointers to points at every component in bodies */
	void setPPointersLoadBodies();						/** Set the pointers to points at every load in bodies */
	//void setPPointersCompJoints();						/** Set the pointers to points at every component in joints */
	//void setPPointersLoadJoints();						/** Set the pointers to points at every load in joints */
	void setPPointersConsJoints();						/** Set the pointers to points at every constraint in joints */	
	void setPPointSurface();							/** Set pointers to surface points in bodies */
	bool checkIfPointProcess(string);				/** Checks if a point belongs to any system within the local process, including root */
	bool checkIfPointSystem(string);				/** Checks if a point belongs to this system */
	bool checkIfPointRootSystem(string);			/** Checks if a point belongs to the root system */
	bool checkIfSystem(string);						/** Checks if a system is in the local process */
	bool checkIfPoint(string);						/** Checks if a point exists on the simulation, based on its system name */
	bool checkIfPointChanged(string);				/** Checks if a point name has changed */
	void addRemoteSharedPoint(string,int);			/** Adds a shared point at a remote process */
	int getProcNum();								/** Get the local process number */
	int getRankNum();								/** Get the local rank */
	void setDofsBodies();							/** Set degrees of freedom at bodies */
	void setDofsJoints();							/** Set degrees of freedom at joints*/
	int countDofs();								/** Counts the total number of degrees of freedom of the system */
	int setTransNumbers(int);						/** Set global numbers for translational dofs */
	int setRotNumbers(int);							/** Set global numbers for rotational dofs */
	int setTRNumbers(int);							/** Set global numbers for translational and rotational dofs */
	void setLocalGlob();							/** Set local-global numbering on components */
	void initializeGlobalDofsBodies(vector<double> &, vector<double> &, double &time);		/** Collect initial dofs from components */
	void initializeGlobalDofsJoints(vector<double> &, vector<double> &, double &time);		/** Collect initial dofs from components */
	void setMatAtComp();							/** Set pointers to material at components */
	Material* getPTMaterial(const string& matName);	/** Return a pointer to the material matName */

	void fillResidual(vector<double>& res_in,
					  vector<double>& q_in,
					  vector<double>& dq_in,					  
					  const double& time, const int&, const double&, const int&, const double&);

	void fillResidual(	vector<double>& residue,
						const vector<double>& q,
						const vector<double>& qdot,
						const vector<double>& qddot,
						const double& time,
						const int&,
						const double&, const int&,
//JC						const int&, const vector<double>&, const vector<double>&,
  						const int&,
						const double&
					);

	void fillResFi(	vector<double>& residue,
					  const vector<double>& q,
					  const vector<double>& qddot,
					  const int&
					  );
	void fillResFq(	vector<double>& residue,
					  const vector<double>& q,
					  const vector<double>& qdot,
					  const int&,
//JC					  const int&, const vector<double>&, const vector<double>&,
   					  const int&,
					  const double&

					  );
	void fillResFt(	vector<double>& residue,
					const double& time,
					const int&,
					const double&, const int&
					  );
//	void fillTangent(sparseMatrix&, vector<double>&, const double& time, const int&);

	void fillTangent(sparseMatrix&, vector<vector<double> >&, vector<vector<double> >&, vector<vector<double> >&, 
					 vector<int>&, vector<int>&, 
					 vector<double>&, vector<double>&,
					 const double& time, const int&);
/*
	void fillTangent(	sparseMatrix& tangent,
						const vector<double>& q,
						const vector<double>& qdot,
						const vector<double>& qddot,
						double partial_qdot,
						double partial_qddot,
						const double& time,
						const int&,
 						const int&,
						const double&
					);
*/
	void fillTangent(sparseMatrix&, vector<vector<double> >&, vector<vector<double> >&, vector<vector<double> >&,
					 vector<int>&, vector<int>&,
					 const vector<double>& q,
					 const vector<double>& qdot,
					 const vector<double>& qddot,
					 double partial_qdot,
					 double partial_qddot,
					 const double& time,
					 const int&,
					 const int&,
					 const double&
					 );
	
	void fillMassMatrix(sparseMatrix& Ma, const vector<double>& q, double time);
	void fillMassMatrix(vector<double>& vMa, const vector<double>& q, double time);
	void fillForce(vector<double>& f,
						const vector<double>& q,
						const vector<double>& qdot,
						const double& time,
						const double&, const int&
					);
	void reportContent();				/** report bodies and systems within the system */
	void reportPointsDofs();			/** report point - dofs */
	void reportPointsDofsFile();		/** report point - dofs */
        void writeFileGraph(ofstream&, string&, string&, string&, bool&, bool&, bool&);		/** write graphics file */
	double getDof(int);
	void setDofsRoot();					/** Set dofs of points at root system */
	vector<string> getPointNames();		/** Get the names of the points that have been already created, from mPoints */
	void updateDofs(vector<double> &, vector<double> &);	/** Update dof vector at point level */
	void setSparsePattern(vector<vector<int> >&);
	void setSparsePatternInertia(vector<vector<int> >&);
	void changePointName(string, string);
	char* getOutGraphFile();
	void initializeCompRes(const bool&, const double&, const int&, const bool);
//	void initializeCompTang(const bool&, const bool&);
	void initializeCompTang(const bool&);
	void setPointSet(string&, string&, vector<double>&);
	map<string,vector<string>,less<string> >* getPointSets()	{return &pointSets;}
	vector<double> getCoordsPoint(string&);	// Get the coordinates of a local point
	void calcEner(double &, double &, double&, double &, double&, double&, double&, int&, vector<double>&, vector<double>&);
	void calcEntropy(double &, double &);
	void calcAM(vector<double> &, vector<double>&, vector<double>&);
	void calcCon(vector<double> &, vector<double>&, vector<double>&, vector<double>&, double&);
	void setIniPosPoint(string&, vector<double>&);
	void setIniPosPointSet(string&, vector<double>&);
	void setIniVelPoint(string&, vector<double>&);
	void setIniVelPointSet(string&, vector<double>&);
	void setIniVelBody(string&, string&, vector<double>&, vector<double>&);
	void setInitialConditions(vector<double>&, vector<double>&);
	void update(const int&, const vector<double>&, const vector<double>&);
	int setInertiaNumbers(int, vector<int>&);
	void setTemperature(double&);
	Body* getPTBodyLocal(string);	/** get a pointer to any body of the local process */
	void writeFileTemp(const double&, ofstream &, string&);
	void setIniPosBody(string&, vector<double>&);
	void setQatComponents(const vector<double>&, const vector<double>&);
	void updateLagMul(void);
	void checkConvLagMul(bool&);
	void calcProjMatrix(sparseMatrix&, double&);
	void fillPhiq(vector<vector<double> >&, int&);
	void setSurfAtBodies();							// Set pointers to surfaces at bodies
//	void setLocalPointsAtBodies();
//	void setBodiesAtJoints();
	void createContacts();
	void addNoDofsPoint(string&, string&, double&);
	void repairChangedPointNames();
	string getNewName(string&);
	void fillMSmooth(sparseMatrix& Ma, const vector<double>& q, double time);
	void fillFSmooth(vector<vector<double> >&, vector<double>&, vector<double>&, double&);
	double getStressValue(int, int);

	int checkComponents(int);
	int setComponents();
	void readTorsionFreeBeam(ifstream&, string&);
    void readRigidBar(ifstream&, string&);

		
	/** Friend functions */
	friend void addSystemName(string& str, const string& strname);
	friend string getStringHead(const string& str);
	friend string getStringTail(const string& str);
	friend void substractSTLVectors(vector<double>& a, vector<double>& b);	 // a-b, a is overwritten	
	friend void addSTLVectors(vector<double>& a, vector<double>& b);	 // a-b, a is overwritten	
	friend void crossProdSTLVectors(vector<double>& a, vector<double>& b);  // a x b, b is overwritten
	friend void printSTLVector(const vector<string>& v);
	
};

#endif
